%% Load entire dataset
% datstruct and IDs from patients with GP struct


clear all; close all; clc;


load ../outputs/dynamics/datastruct_v1v2v5_it960.mat
ids = ids(1:960); % due to an error in the script, the last 8 patients' results were not saved.

run('../Toolboxes/GP_v3/startup.m'); % Add directories of GP_v3 toolbox to path

% Plot example patient...
DOPLOT = 1;
pt = 800;   % example patient index
vt = 5;     % vital sign [HR, PP, LA, UO, MAP]
if DOPLOT
    figure;
    zz = datastruct{ids(pt)}.gp{vt}.xtest; 
    mu = datastruct{ids(pt)}.gp{vt}.ytest(:,1); 
    s2 = datastruct{ids(pt)}.gp{vt}.ytest(:,2);
    x = datastruct{ids(pt)}.gp{vt}.xtrain;
    y = datastruct{ids(pt)}.gp{vt}.ytrain;
    f = [mu+2*sqrt(s2); flipdim(mu-2*sqrt(s2),1)];
    fill([zz; flipdim(zz,1)], f, [7 7 7]/8); alpha(0.25); hold on;
    plot(zz, mu,'--','Color','k'); plot(x,y,'*','Color','r');
    hold off; xlabel('x, input'); ylabel('y, output');
    title('GP-Regression');
end

% Check if data are all there
for i = 1 : length(ids)
    if ~isfield(datastruct{ids(i)},'gp')
        fprintf('\nPROBLEM!\n');
    end
end

%% Perform clustering with Bhat Distance...
% Pick variables to perform joint clustering

% variab = [1]; % Perform comparison using these variables
% variabname = '';
% for v = variab
%     variabname = sprintf('%sv%d',variabname,v);
% end
% MATRIX = -1*ones(length(ids));
% for i = 1 : length(ids)
%     clear Y1;
%     for v = 1 : length(variab)
%         Y1{v} = datastruct{ids(i)}.gp{variab(v)}.ytest;
%     end
%     for j = i : length(ids)
%         % extract time series
%         clear Y2;
%         for v = 1 : length(variab)
%             Y2{v} = datastruct{ids(j)}.gp{variab(v)}.ytest;
%         end
%         MATRIX(i,j) = computegpdistance(Y1,Y2);
%         MATRIX(j,i) = MATRIX(i,j);
%     end
%     fprintf('\nComplete %.1f',i/length(ids)*100);
% end
% fprintf('\n');
% disp('');
% save(sprintf('../outputs/dynamics/matrix_%s.mat',variabname),'MATRIX');
% 
% MATRIX1 = MATRIX*2 - 1;
% clustergram(MATRIX1);
% 
% Z1 = convert_matrix2pdist(MATRIX);
% Z = linkage(Z1);

%% Perform clustering with inferred hyperparameters
% Get hyperparameters and perform clustering with k-means

variab = [1 2 5]; % Perform comparison using these variables
featureSpace = zeros(length(ids),8*length(variab));
% FeatureList = [cov4 mean std mean_ mean+]
mort = -1*ones(length(ids),1);
for i = 1 : length(ids)
    n = 1;
    for v = 1 : length(variab)
        featureSpace(i,n:n+3) = datastruct{ids(i)}.gp{variab(v)}.hyp_post.cov;
        featureSpace(i,n+4) = mean(datastruct{ids(i)}.gp{variab(v)}.ytest(:,1));
        featureSpace(i,n+5) = std(datastruct{ids(i)}.gp{variab(v)}.ytest(:,1));
        N = length(datastruct{ids(i)}.gp{variab(v)}.ytest(:,1)); N_2 = round(N/2);
        featureSpace(i,n+6) = mean(datastruct{ids(i)}.gp{variab(v)}.ytest(1:N_2,1));
        featureSpace(i,n+7) = mean(datastruct{ids(i)}.gp{variab(v)}.ytest(N_2:N,1));
        n = n + 8;
    end
    mort(i,1) = datastruct{ids(i)}.mort;
    fprintf('\nComplete %.1f',i/length(ids)*100);
end
% remove survivor patients...
% ind = find(mort == 0);
fprintf('\n');
featureSelect = [1 3 7 8];
featureSelect = [featureSelect featureSelect+8 featureSelect+16];
Xa = featureSpace(:,featureSelect);
X = Xa-repmat(mean(Xa),size(Xa,1),1);
X = X./repmat(std(Xa),size(Xa,1),1);
% X = Xa;
k = [5 10 20];
cdata = [.3 .2 .5];
figure;
for i = 1 : length(k)
    [IDX{i},C{i}] = kmeans(X,k(i));
    [count,xcentres] = hist(IDX{i},1:k(i));
    count = count/sum(count);
    countD = zeros(size(count));
    for c = 1 : k(i)
        ind = find(IDX{i} == c);
        countD(c) = sum(mort(ind))/length(ind);
    end
    subplot(2,3,i);
    bar(xcentres,count,'FaceColor',cdata,'EdgeColor','k');
    ylabel('Mixture Proportion'); title('Mixture Proportion per Cluster');
    xlim([.5 k(i)+.5]);
    subplot(2,3,i+3);
    bar(xcentres,countD,'FaceColor',cdata,'EdgeColor','k');
    hold on; plot([0 k(i)+1], [sum(mort)/length(ids) sum(mort)/length(ids)], 'k',...
        'LineWidth',2);
    ylabel('Mortality Proportion'); title('Mortality proportion per Cluster');
    xlim([.5 k(i)+.5]); ylim([0 .8])
end

%% Perform clustering analysis...

% For survivors group
c = 13; % 
ind = find(IDX{3} == c);
ind =ind(:)';

for pt = ind
    figure;
    subplot(2,2,1)
    bar(xcentres,countD,'FaceColor',cdata,'EdgeColor','k');
    hold on; plot([0 k(i)+1], [sum(mort)/length(ids) sum(mort)/length(ids)], 'k',...
        'LineWidth',2);
    ylabel('Mortality Proportion'); title('Mortality proportion per Cluster');
    xlim([.5 k(i)+.5]); ylim([0 .8])
    for i = 1 : 3
        subplot(2,2,i+1); vt = variab(i);
        zz = datastruct{ids(pt)}.gp{vt}.xtest; 
        mu = datastruct{ids(pt)}.gp{vt}.ytest(:,1); 
        s2 = datastruct{ids(pt)}.gp{vt}.ytest(:,2);
        x = datastruct{ids(pt)}.gp{vt}.xtrain;
        y = datastruct{ids(pt)}.gp{vt}.ytrain;
        f = [mu+2*sqrt(s2); flipdim(mu-2*sqrt(s2),1)];
        fill([zz; flipdim(zz,1)], f, [8 5 5]/8); alpha(0.25); hold on;
        plot(zz, mu,'--','Color','k'); plot(x,y,'*','Color','r');
        hold off; xlabel('x, input'); ylabel('y, output');
        title(sprintf('Death? %d',datastruct{ids(pt)}.mort));
    end
    pause;
end